Biologically informed NeuralODEs for genome-wide regulatory dynamics

Background Gene regulatory network (GRN) models that are formulated as ordinary differential equations (ODEs) can accurately explain temporal gene expression patterns and promise to yield new insights into important cellular processes, disease progression, and intervention design. Learning such gene regulatory ODEs is challenging, since we want to predict the evolution of gene expression in a way that accurately encodes the underlying GRN governing the dynamics and the nonlinear functional relationships between genes. Most widely used ODE estimation methods either impose too many parametric restrictions or are not guided by meaningful biological insights, both of which impede either scalability, explainability, or both. Results We developed PHOENIX, a modeling framework based on neural ordinary differential equations (NeuralODEs) and Hill-Langmuir kinetics, that overcomes limitations of other methods by flexibly incorporating prior domain knowledge and biological constraints to promote sparse, biologically interpretable representations of GRN ODEs. We tested the accuracy of PHOENIX in a series of in silico experiments, benchmarking it against several currently used tools. We demonstrated PHOENIX’s flexibility by modeling regulation of oscillating expression profiles obtained from synchronized yeast cells. We also assessed the scalability of PHOENIX by modeling genome-scale GRNs for breast cancer samples ordered in pseudotime and for B cells treated with Rituximab. Conclusions PHOENIX uses a combination of user-defined prior knowledge and functional forms from systems biology to encode biological “first principles” as soft constraints on the GRN allowing us to predict subsequent gene expression patterns in a biologically explainable manner. Supplementary Information The online version contains supplementary material available at 10.1186/s13059-024-03264-0.

Additional File 1: Figure S2 We measured the marginal contribution of PHOENIX's architecture and incorporation of prior information by comparing against the baseline contributions of out-of-the-box NeuralODE models with three different activation functions, across both in silico dynamical systems SIM350 and SIM690, after subjecting the data to different amounts of noise.We display the performance of the out-of-the-box models, as well as PHOENIX (λ prior tuned using the validation set) and its unregularized version (λ prior = 0), in terms of how well held out time points from noiseless test set trajectories could be predicted after training on trajectories from different noise settings.The entire procedure was repeated five times to generate average mean-squared error (MSE) values and error bars.
Additional File 1: Table S1 Assessing the contribution of priors to PHOENIX explainability, as well as the effect of prior 1 For in silico experiments (SIM350, SIM690) under 5% noise, the prior knowledge model P * was rendered "incomplete" by setting a certain percentage of nonzero prior interactions to 0 (see Additional File 2: Section 4.2).
2 For each scenario, an AUC was calculated by aligning the entire ground truth GRN against a network extracted from a PHOENIX model trained using the incomplete prior in question.AUC values were also separately calculated for subsets of the ground truth GRN edges based on whether or not the edge was withheld from P * .
Additional File 1: Figure S3 We extracted encoded GRNs from the trained PHOENIX models and the best-performing out-of-the-box NeuralODE models, for both in silico dynamical systems SIM350 (A) and SIM690 (B) across all noise settings.Here we display noise settings 40%, 80% and 100%.These results are provided in addition to the noise settings of 0%, 5%, 10%, and 20% that have been shown in Figure 3 of the main paper.We compared these GRN estimates to the corresponding ground truth GRNs used to formulate SIM350 and SIM690, and obtained AUC values as well as out-degree correlations (ρout).We also reverse-engineered a metric (Cmax) to inform how sparsely PHOENIX had inferred the dynamics (see Additional File 2: Section 2).Furthermore, we used these Cmax values to obtain optimal true positive and true negative rates (TPRmax and TNRmax) that were independent of any cutoff value, allowing us to compare between "best possible" networks across all settings.
Additional File 1: Figure S4 We investigated how PHOENIX used prior information to sparsify the learned dynamics by computing the encoded 350 × 350 dynamics matrices D (see Additional File 2: Section 2), across multiple noise settings in SIM350.We display D for both PHOENIX as well as its unregularized (λ prior = 0) version as heatmaps of effect size.Unlike its unregularized counterpart, PHOENIX identifies core elements of the dynamics even at relatively high noise levels.
Additional File 1: Figure S5 PHOENIX uses gene-specific multipliers υ ∈ R n (see Model Formulation in main paper) to simplify the representation of steady state for genes without upstream transcription factors = 0, ∀t .In order to explore the contribution of these gene-specific multipliers towards the PHOENIX model, we tested PHOENIX with (A) and without (B) the multipliers using gene expression data originating from the in silico dynamical system SIM690.We used SIM690 to simulate the temporal gene expression patterns for its 690 genes.Each simulated trajectory consisted of five time points (t = 0, 2, 3, 7, 9) and was subjected to varying levels of Gaussian noise ( noise σ mean = 0%, 5%, 10%, and 20% ).We also corrupted the user-defined prior models up to an amount commensurate with the noise level.We display across all noise settings both observed and predicted trajectories for four arbitrary genes that had true flat trajectories i.e dg i (t) dt = 0, ∀t in SIM690.
Additional File 1: Table S3 In order to incorporate prior domain knowledge, PHOENIX uses an adjacency matrix A of likely network structure based on user-defined biological insights (including experimentally validated interactions, motif map of promoter targets, etc.).Here, A ij ∈ {+1, −1, 0} representing an activating, repressive, or no prior interaction, respectively.But for real world scenario, the signs (activating/repressive) of prior interactions are often unknown.So we wanted to check whether formulating A simply based on prior interaction existence, A ij ∈ {+1, 0}, would suffice.Below we show the results of testing the importance of providing the correct prior signs to PHOENIX using our in silico setup where the ground truth signs were always known.We tested multiple forms of regularized PHOENIX models and report results both in terms predictive accuracy (MSE on test set) and explainability (AUC from comparing inferred dynamics to ground truth GRNs).
2 PHX + = Regularized PHOENIX where λ prior was chosen based on the validation set, but both activating and repressive edges in the prior were set to +1. "No prior interaction" was kept to 0. This emulated a setting where we didn't have prior knowledge of these signs (either activating vs repressive prior edge).
3 PHX = Regularized PHOENIX where λ prior was chosen based on the validation set.Each entry in the prior was correctly set to +1 or -1, depending on whether the prior edge was activating or repressive."No prior interaction" was set to 0.
Additional File 1: 1 PHX = Regularized PHOENIX where λ prior was chosen based on the validation set.
3 OOTB = Out-of-the-box NeuralODE, resembling how plain NeuralODEs are typically used for this problem.We note that the method RNAForecaster uses an OOTB approach with ReLU activation, but the results here are for the tanh activation function that had better performance than both sigmoid and ReLU on the validation set.
4 PRESC = PRESCIENT, where we used the validation set to optimize the value of k dim , which is the number of neurons use in PRESCIENT's hidden layer.
5 For Dynamo we used the validation set to optimize the sparsity regularization penalty (λ), as well as the number of kernel basis functions use to approximate the vector field (M ) 6 RNODE = RNA-ODE, where we used the validation set to optimize the number of trees used in the random forest function 7 In DeepVelo the encoder and the decoder consisted of four dense layers (size 4p for the intermediate layers and size p for the latent layer) with ReLU activation.We used the validation set to optimize p.
Additional File 1: PBA ostensibly ignores oscillatory gene expression dynamics (the cell cycle), and hence may not be flexible enough.Also, we already benchmarked against another method (PRESCIENT ), that assumes dynamics to be the gradient of a scalar potential.

PathReg
Uses an out-of-the-box NeuralODE with an Exponential Linear Unit (ELU) activation function and enforces both weight and feature sparsity at the path-level via a differentiable L 0 -based regularizer in the loss function.A matrix of non-negative stochastic gates regularizes the probability of any path throughout the entire network contributing to a given output.This constrains the number of inputoutput paths.
PathReg was primarily developed to induce sparsity in NeuralODEs, and is similar in principle to the out-of-the-box models that we already benchmarked.Consequently, is has only been demonstrated to scale up to 529 highly variable dimensions/genes.

PHOENIX (Supplemental Results
Additional File 1: Table S8  Additional File 1: Figure S6 We applied PHOENIX (λ prior = 0.05) to 2 technical replicates of gene expression of 3551 genes each, collected across 24 time points in a yeast cell-cycle time course.We trained on 40 transition pairs, used 3 for validation, and tested predictive accuracy on the remaining 3. We display both observed and predicted trajectories for ACT1 and TDH1, where the predicted trajectories are extrapolations into future time points based on just initial values (gene expression at t = 0).
Additional File 1: Table S9 PHOENIX uses gene-specific multipliers υ ∈ R n (see Model Formulation in main paper) to simplify the representation of steady state for genes without upstream transcription factors dg i (t) dt = 0, ∀t .Below we report the number of genes for which υ i ≤ 0 (and hence gene i has a flat trajectory) for the PHOENIX model fitted on the yeast cell cycle data.

Dataset Number of genes modelled Number of genes
Additional File 1: Figure S10 For the breast cancer data, we depict the R 2 performance (y-axis) of a PHOENIX model trained on all genes, when considering only the N most variable genes (x-axis) for evaluation.We further evaluate the impact of predicting the entire trajectory based on the initial gene expression value (t 0 , blue) versus predicting each expression value based on its immediately previous time point (red).

PHOENIX (Supplemental Results
Additional File 1: Table S10 Detailed results of permutation tests (see Methods 3.4) used to obtain pathway influence scores from trained PHOENIX models and the Reactome pathway database.The mean (µ 0 ) and standard deviation (σ 0 ) of each permutation test null distribution, over K = 1000 permutations, are tabulated for each subset (ng) in question, and the corresponding z-scores are visualized in Figure 5. Missing values indicate that the all genes involved in that pathway were excluded from that particular subset based on low variability in expression.
Reactome accession IDs for all the pathways below have been provided in Additional File 1: Table S11 for reproducibility.
Reactome pathway ng = 500 genes ng = 2000 genes ng = 4000 genes ng = 11165 genes 1 The full data set consisted of ng = 11165 genes.We also fit PHOENIX to smaller subsets of genes ng = 500, 2000, 4000, where we subsetted the full data set to only the ng most variable genes in the pseudotrajectory.
3 For each ng, we computed inferred influence scores for individual genes based on perturbation analyses on the fitted PHOENIX model (see Methods 3.3).We also computed harmonic centralities of the genes based on a validation network describing ChIP-binding data.We then did a binary assignment, where the top 10% (10% × ng) genes with highest inferred influence were labelled "predicted influential" and the remaining (90%×ng) were "predicted non-influential."We measured concordance as the sensitivity T P T P +F N with which these predicted labels recovered the "truly influential" genes (those genes with non-zero harmonic centrality in the ChIP validation network).
Additional File 1: Table S14 We compared the performance of PHOENIX to other methods of regulatory dynamics estimation on a pseudotrajectory of 186 breast cancer samples (ordered along subsequent "pseudotimepoints") consisting of ng = 11165 genes.The data was processed for model fitting via the steps described in Methods 3.2.Snapshot-based methods (Dynamo, RNA-ODE, Deepvelo) require RNA velocity at every time point as an additional input.Given that this information was not available in the data set, we estimated RNA velocity using a method of finite differences applied to smooth splines through the expression trajectories (see Additional File 2: Section 3.1).Once each model was trained, we measured predictive accuracy in a step-wise fashion; for N ranging from 50 to 11165, we calculated the test set MSEs (as described in Methods 3.2) of only the top-N most varying genes in the dataset.Additional File 1: Table S16 Training results from PHOENIX models applied to the B-cell RNASeq datasets.We fitted two seperate PHOENIX models, one for each condition of untreated and treated with Rituximab.Furthermore, PHOENIX uses gene-specific multipliers υ ∈ R n (see Model Formulation in main paper) to simplify the representation of steady state for genes without upstream transcription factors dg i (t) dt = 0, ∀t ; so we report the number of genes for which υ i ≤ 0 (and hence gene i has a flat trajectory) for each of two models.

Condition
misspecification, in both in silico experiments and real experimental data Additional File 1: TableS2Assessing the effect on explainability of withholding complete subsets of information from the prior (P * ), in in silico experiments.Gaussain noise of 5% was added to all expression values.
Table S4 Benchmarking PHOENIX against other methods in silico, in terms of trajectory-recovery performance (performance metric is MSE on test set) TableS5Benchmarking PHOENIX against other methods in silico, in terms of explainability.A GRN describing the inferred dynamics was extracted for each fitted model, and was compared to the ground truth GRN to calculate an AUC DeepVelo provides a function for estimating gene correlation networks from simulated retrograde trajectories.We generated n = 200 retrograde trajectories using 200 randomly generated initial conditions (see Additional File 2: Section 3.3)Additional File 1: TableS6Benchmarking PHOENIX against other methods in silico, in terms of sparsity of inferred dynamics.Sparsity was calculated here as the average out-degree of the extracted GRN from each of the inferred dynamical systems Additional File 1: TableS7Details about additional black-box methods for estimating gene expression dynamics that we excluded from our in silico benchmarking experiments.We provide details about each method, and some reasoning behind its exclusion 1The architecture of PHOENIX allowed GRNs to be extracted from both regularized and unregularized versions using a very simple and efficient algorithm (see Additional File 2: Section 2) 2 GRNs were extracted from out-of-the-box NeuralODEs, using sensitivity analyses (see Additional File 2: Section 3.2) 3 PRESCIENT does not provide a straightforward means for extracting a GRN describing the dynamics 4 Dynamo provides a function that calculates the Jacobian matrix of the fitted model at any data point.We used this function to calculate an average Jacobian matrix across several randomly generated data points and estimate a GRN (see Additional File 2: Section 3.3) 5 RNA-ODE provides a function for GRN inference by ranking regulatory links using estimated effect sizes (see Additional File 2: Section 3.3) 6 Qualitative comparison of black-box methods for estimating gene expression dynamics that we benchmarked against PHOENIX in silico.We provide details about inputs, learning algorithms, and key performance metrics Additional File 1: TableS12PHOENIX uses gene-specific multipliers υ ∈ R n (see Model Formulation in main paper) to simplify the representation of steady state for genes without upstream transcription factors dg i (t) dt = 0, ∀t .Below we report the number of genes for which υ i ≤ 0 (and hence gene i has a flat trajectory) for each of the PHOENIX models fitted on the breast cancer data.
Additional File 1: TableS15Top 50 most changing regulators discovered by PHOENIX in B cells after Rituximab treatment as measured by log-fold change.
#Genes modelled val.set MSE train set R 2 #Genes with υ i ≤ 0